pacman::p_load(maptools, sf, raster, spatstat, tmap, tidyverse, funModeling)Take-Home Exercise 1
1 Overview
1.1 Context
Water is an essential resource for people. The health of people depends on having access to clean water. It ensures peace and security, a clean environment, a sustainable economy, and a reduction in poverty. Yet more than 40% of the world’s population lacks access to enough clean water. According to UN-Water, 1.8 billion people will be residing in nations or areas with a severe water shortage by 2025. Water scarcity is a serious threat to several sectors, including food security. Around 70% of the world’s freshwater resources are used for agriculture.
1.2 Objectives
Discover the geographical distribution of functional and non-function water points and their co-locations if any in Osun State, Nigeria.
Exploratory Spatial Data Analysis (ESDA)
Derive kernel density maps of functional and non-functional water points. Using appropriate tmap functions
Display the kernel density maps on openstreetmap of Osub State, Nigeria.
Describe the spatial patterns revealed by the kernel density maps. Highlight the advantage of kernel density map over point map.
Second-order Spatial Point Patterns Analysis
With reference to the spatial point patterns observed in ESDA:
Formulate the null hypothesis and alternative hypothesis and select the confidence level.
Perform the test by using appropriate Second order spatial point patterns analysis technique.
With reference to the analysis results, draw statistical conclusions.
Spatial Correlation Analysis
Investigate if the spatial distribution of functional and non-functional water points are independent from each other.
Formulate the null hypothesis and alternative hypothesis and select the confidence level.
Perform the test by using appropriate Second order spatial point patterns analysis technique. With reference to the analysis results, draw statistical conclusions.
1.3 Data
| Type | Name | Format | Description | Source |
|---|---|---|---|---|
| Aspatial | WPdx+ | csv | Locations of water points | WPdx Global Data Repositories |
| Geospatial | geoBoundaries | SHP | geoBoundaries data of Nigeria | geoBoundaries |
2 Importing and loading packages
For this exercise, we’ll be using the following packages:
sf: Manage and process vector-based geospatial data in R
spatstat: Perform 1st & 2nd spatial point patterns analysis + kernel density
raster: Reads, writes, manipulates, analyses, model of grid spatial data convert image output generate by spatstat into raster format
maptools: Convert Spatial objects into ppp format of spatstat
tmap: Plotting cartographic quality static point patterns maps or interactive maps by using leaflet API
tidyverse: a collection of functions, data, and documentation that extends the capabilities of base R.
funModeling: Plot charts for easier interpretations
We will be using p_load function of pacman package to install and load required packages.
3 Spatial Data Wrangling
3.1 Data Geospatial
We will be using st_read function to read our geospatial data.
NGA <- st_read(dsn = "data/geospatial/",
layer = "nga_admbnda_adm2_osgof_20190417")%>%
st_transform(crs = 26392)Reading layer `nga_admbnda_adm2_osgof_20190417' from data source
`C:\kt-x\is415-GAA\Take-Home_Ex\Take-Home_Ex01\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 774 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2.668534 ymin: 4.273007 xmax: 14.67882 ymax: 13.89442
Geodetic CRS: WGS 84Lets start by selecting specific fields that will be helpful in this exercise and select Osun state using the filter() function.
NGA <- NGA %>%
select(c(3:4, 8:9)) %>%
filter(ADM1_EN == "Osun")glimpse(NGA)Rows: 30
Columns: 5
$ ADM2_EN <chr> "Aiyedade", "Aiyedire", "Atakumosa East", "Atakumosa West",…
$ ADM2_PCODE <chr> "NG030001", "NG030002", "NG030003", "NG030004", "NG030005",…
$ ADM1_EN <chr> "Osun", "Osun", "Osun", "Osun", "Osun", "Osun", "Osun", "Os…
$ ADM1_PCODE <chr> "NG030", "NG030", "NG030", "NG030", "NG030", "NG030", "NG03…
$ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((213526.6 34..., MULTIPOLYGON (…3.2 Importing WPdx
We will be using read_csv() function to read the WPdx file and using filter function to only show Nigeria country data and only for Osun state.
wp_nga <- read_csv("data/aspatial/WPdx.csv") %>%
filter(`#clean_country_name` == "Nigeria") %>%
filter(`#clean_adm1` == "Osun")3.2.1 Convert water point data into sf point features
wp_nga$Geometry = st_as_sfc(wp_nga$`New Georeferenced Column`)
wp_nga# A tibble: 5,557 × 71
row_id `#source` #lat_…¹ #lon_…² #repo…³ #stat…⁴ #wate…⁵ #wate…⁶ #wate…⁷
<dbl> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr>
1 429123 GRID3 8.02 5.06 08/29/… Unknown <NA> <NA> Tapsta…
2 70566 Federal Minis… 7.32 4.79 05/11/… No Protec… Well Mechan…
3 70578 Federal Minis… 7.76 4.56 05/11/… No Boreho… Well Mechan…
4 66401 Federal Minis… 8.03 4.64 04/30/… No Boreho… Well Mechan…
5 422190 GRID3 7.87 4.88 08/29/… Unknown <NA> <NA> Tapsta…
6 422064 GRID3 7.7 4.89 08/29/… Unknown <NA> <NA> Tapsta…
7 65607 Federal Minis… 7.89 4.71 05/12/… No Boreho… Well Mechan…
8 68989 Federal Minis… 7.51 4.27 05/07/… No Boreho… Well <NA>
9 67708 Federal Minis… 7.48 4.35 04/29/… Yes Boreho… Well Mechan…
10 66419 Federal Minis… 7.63 4.50 05/08/… Yes Boreho… Well Hand P…
# … with 5,547 more rows, 62 more variables: `#water_tech_category` <chr>,
# `#facility_type` <chr>, `#clean_country_name` <chr>, `#clean_adm1` <chr>,
# `#clean_adm2` <chr>, `#clean_adm3` <chr>, `#clean_adm4` <chr>,
# `#install_year` <dbl>, `#installer` <chr>, `#rehab_year` <lgl>,
# `#rehabilitator` <lgl>, `#management_clean` <chr>, `#status_clean` <chr>,
# `#pay` <chr>, `#fecal_coliform_presence` <chr>,
# `#fecal_coliform_value` <dbl>, `#subjective_quality` <chr>, …3.2.2 Convert this df into sf object
wp_sf <- st_sf(wp_nga, crs=4326)
wp_sfSimple feature collection with 5557 features and 70 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 4.032004 ymin: 7.060309 xmax: 5.06 ymax: 8.061898
Geodetic CRS: WGS 84
# A tibble: 5,557 × 71
row_id `#source` #lat_…¹ #lon_…² #repo…³ #stat…⁴ #wate…⁵ #wate…⁶ #wate…⁷
* <dbl> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr>
1 429123 GRID3 8.02 5.06 08/29/… Unknown <NA> <NA> Tapsta…
2 70566 Federal Minis… 7.32 4.79 05/11/… No Protec… Well Mechan…
3 70578 Federal Minis… 7.76 4.56 05/11/… No Boreho… Well Mechan…
4 66401 Federal Minis… 8.03 4.64 04/30/… No Boreho… Well Mechan…
5 422190 GRID3 7.87 4.88 08/29/… Unknown <NA> <NA> Tapsta…
6 422064 GRID3 7.7 4.89 08/29/… Unknown <NA> <NA> Tapsta…
7 65607 Federal Minis… 7.89 4.71 05/12/… No Boreho… Well Mechan…
8 68989 Federal Minis… 7.51 4.27 05/07/… No Boreho… Well <NA>
9 67708 Federal Minis… 7.48 4.35 04/29/… Yes Boreho… Well Mechan…
10 66419 Federal Minis… 7.63 4.50 05/08/… Yes Boreho… Well Hand P…
# … with 5,547 more rows, 62 more variables: `#water_tech_category` <chr>,
# `#facility_type` <chr>, `#clean_country_name` <chr>, `#clean_adm1` <chr>,
# `#clean_adm2` <chr>, `#clean_adm3` <chr>, `#clean_adm4` <chr>,
# `#install_year` <dbl>, `#installer` <chr>, `#rehab_year` <lgl>,
# `#rehabilitator` <lgl>, `#management_clean` <chr>, `#status_clean` <chr>,
# `#pay` <chr>, `#fecal_coliform_presence` <chr>,
# `#fecal_coliform_value` <dbl>, `#subjective_quality` <chr>, …3.2.3 Transforming into Nigeria projected coordinate system
wp_sf <- wp_sf %>%
st_transform(crs = 26392)3.2.4 EDA for Water Point Data
Next, we want to first have a visual of the water point data by creating a frequency chart. (Reference to In-class exercise 2)
freq(data = wp_sf,
input = '#status_clean')
#status_clean frequency percentage cumulative_perc
1 Functional 2319 41.73 41.73
2 Non-Functional 2008 36.13 77.86
3 <NA> 748 13.46 91.32
4 Functional but needs repair 248 4.46 95.78
5 Non-Functional due to dry season 151 2.72 98.50
6 Functional but not in use 63 1.13 99.63
7 Abandoned 15 0.27 99.90
8 Abandoned/Decommissioned 5 0.09 100.00Seems like there are several different categories of water points. With reference to Objectives, we will be focusing on Functional and Non-Functional while still taking into account of the unknown water points. But first, let’s rename this spatial frame.
wp_sf_nga <- wp_sf %>%
rename(status_clean = '#status_clean') %>%
select(status_clean) %>%
mutate(status_clean = replace_na(
status_clean, "unknown"))After that, we filter and assign the wp_functional
wp_functional <- wp_sf_nga %>%
filter(status_clean %in%
c("Functional",
"Functional but not in use",
"Functional but needs repair"))
freq(data = wp_functional,
input = 'status_clean')
status_clean frequency percentage cumulative_perc
1 Functional 2319 88.17 88.17
2 Functional but needs repair 248 9.43 97.60
3 Functional but not in use 63 2.40 100.00Sames goes for the Non-Functional
wp_nonfunctional <- wp_sf_nga %>%
filter(status_clean %in%
c("Abandoned/Decommissioned",
"Abandoned",
"Non-Functional due to dry season",
"Non-Functional",
"Non functional due to dry season"))
freq(data = wp_nonfunctional,
input = 'status_clean')
status_clean frequency percentage cumulative_perc
1 Non-Functional 2008 92.15 92.15
2 Non-Functional due to dry season 151 6.93 99.08
3 Abandoned 15 0.69 99.77
4 Abandoned/Decommissioned 5 0.23 100.00Sames goes for the Unknown
wp_unknown <- wp_sf_nga %>%
filter(status_clean == "unknown")
freq(data = wp_unknown,
input = 'status_clean')
status_clean frequency percentage cumulative_perc
1 unknown 748 100 1003.2.5 Performing Point-in-Polygon Count
Next, we want to find out the number of total, functional, nonfunctional and unknown water points in each LGA. First, it identifies the functional water points in each LGA by using st_intersects() of sf package. Next, length() is used to calculate the number of functional water points that fall inside each LGA.
NGA_wp <- NGA %>%
mutate(`total_wp` = lengths(
st_intersects(NGA, wp_sf_nga))) %>%
mutate(`wp_functional` = lengths(
st_intersects(NGA, wp_functional))) %>%
mutate(`wp_nonfunctional` = lengths(
st_intersects(NGA, wp_nonfunctional))) %>%
mutate(`wp_unknown` = lengths(
st_intersects(NGA, wp_unknown)))4 Geospatial Data wrangling
Moving to data wrangling on geospatial data, first, we need to convert simple feature data frame to sp spatial class.
4.1 Converting sf data frames to sp’s Spatial* class
Convert to sp object/class using as_Spatial(). Take a look at the following 2 tabs (wp_sc, NGA_sc), notice the properties is “SpatialPolygonsDataFrame”.
#wp_sc <- as_Spatial(wp_sf)
wp_sc_functional <- as_Spatial(wp_functional)
wp_sc_nonfunctional <- as_Spatial(wp_nonfunctional)
NGA_sc <- as_Spatial(NGA)#wp_sc
wp_sc_functionalclass : SpatialPointsDataFrame
features : 2630
extent : 177285.9, 290751, 343128.1, 450859.7 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs
variables : 1
names : status_clean
min values : Functional
max values : Functional but not in use wp_sc_nonfunctionalclass : SpatialPointsDataFrame
features : 2179
extent : 180539, 290616, 340054.1, 450780.1 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs
variables : 1
names : status_clean
min values : Abandoned
max values : Non-Functional due to dry season NGA_scclass : SpatialPolygonsDataFrame
features : 30
extent : 176503.2, 291043.8, 331434.7, 454520.1 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs
variables : 4
names : ADM2_EN, ADM2_PCODE, ADM1_EN, ADM1_PCODE
min values : Aiyedade, NG030001, Osun, NG030
max values : Osogbo, NG030030, Osun, NG030 4.2 Converting the Spatial* class into generic sp format
WHY? Because spatstat requires the analytical data in ppp object form. There is no direct way to convert a Spatial* classes into ppp object. We need to convert the Spatial classes* into Spatial object first.
Starting with Converting the Spatial* classes into generic sp objects. Take a look at the following 2 tabs (wp_sc, NGA_sc), notice the properties is “SpatialPoints” and “SpatialPolygons” respectively.
#wp_sp <- as(wp_sc, "SpatialPoints")
wp_sp_functional <- as(wp_sc_functional, "SpatialPoints")
wp_sp_nonfunctional <- as(wp_sc_nonfunctional, "SpatialPoints")
nga_sp <- as(NGA_sc, "SpatialPolygons")#wp_sp
wp_sp_functionalclass : SpatialPoints
features : 2630
extent : 177285.9, 290751, 343128.1, 450859.7 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs wp_sp_nonfunctionalclass : SpatialPoints
features : 2179
extent : 180539, 290616, 340054.1, 450780.1 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs nga_spclass : SpatialPolygons
features : 30
extent : 176503.2, 291043.8, 331434.7, 454520.1 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 4.3 Converting the generic sp format into spatstat’s ppp format
Next, we will use as.ppp() function of spatstat to convert the spatial data into spatstat’s ppp object format.
#wp_ppp <- as(wp_sp, "ppp")
wp_ppp_functional <- as(wp_sp_functional, "ppp")
wp_ppp_nonfunctional <- as(wp_sp_nonfunctional, "ppp")
#wp_ppp
wp_ppp_functionalPlanar point pattern: 2630 points
window: rectangle = [177285.9, 290750.96] x [343128.1, 450859.7] unitswp_ppp_nonfunctionalPlanar point pattern: 2179 points
window: rectangle = [180538.96, 290616] x [340054.1, 450780.1] units#plot(wp_ppp)
par(mfrow=c(1,2))
plot(wp_ppp_functional)
plot(wp_ppp_nonfunctional)
4.3.1 Check for duplicates
#any(duplicated(wp_ppp))
any(duplicated(wp_ppp_functional))[1] FALSEany(duplicated(wp_ppp_nonfunctional))[1] FALSEIf so, check for the number of duplicates
sum(multiplicity(wp_ppp_functional) > 1)[1] 0sum(multiplicity(wp_ppp_nonfunctional) > 1)[1] 0To resolve this problem, we will be using the jittering approach, which will add a small perturbation to the duplicate points so that they do not occupy the exact same space.
#wp_ppp_jit <- rjitter(wp_ppp,
# retry=TRUE,
# nsim=1,
# drop=TRUE)
wp_ppp_functional_jit <- rjitter(wp_ppp_functional,
retry=TRUE,
nsim=1,
drop=TRUE)
wp_ppp_nonfunctional_jit <- rjitter(wp_ppp_nonfunctional,
retry=TRUE,
nsim=1,
drop=TRUE)Then, check again for duplicates
any(duplicated(wp_ppp_functional_jit))[1] FALSEany(duplicated(wp_ppp_nonfunctional_jit))[1] FALSE4.4 Create owin object
owin object is designed to represent this polygonal region. We will be using to convert nigeria Spatial Polygon object into owin object of spatstat.
nga_owin <- as(nga_sp, "owin")plot(nga_owin)
summary(nga_owin)Window: polygonal boundary
30 separate polygons (no holes)
vertices area relative.area
polygon 1 204 766084000 0.08870
polygon 2 81 304399000 0.03520
polygon 3 97 465688000 0.05390
polygon 4 124 373051000 0.04320
polygon 5 60 149473000 0.01730
polygon 6 84 144820000 0.01680
polygon 7 50 102243000 0.01180
polygon 8 72 216002000 0.02500
polygon 9 112 269897000 0.03130
polygon 10 125 365142000 0.04230
polygon 11 83 111191000 0.01290
polygon 12 126 192557000 0.02230
polygon 13 219 904397000 0.10500
polygon 14 174 741131000 0.08580
polygon 15 81 138742000 0.01610
polygon 16 65 119452000 0.01380
polygon 17 90 280205000 0.03240
polygon 18 69 69814600 0.00808
polygon 19 69 42727500 0.00495
polygon 20 49 30458800 0.00353
polygon 21 62 263505000 0.03050
polygon 22 93 438930000 0.05080
polygon 23 87 274127000 0.03170
polygon 24 105 509979000 0.05910
polygon 25 98 292058000 0.03380
polygon 26 64 327765000 0.03800
polygon 27 133 108945000 0.01260
polygon 28 122 462169000 0.05350
polygon 29 94 109715000 0.01270
polygon 30 95 61239800 0.00709
enclosing rectangle: [176503.22, 291043.82] x [331434.7, 454520.1] units
(114500 x 123100 units)
Window area = 8635910000 square units
Fraction of frame area: 0.6134.4.1 Combining point events object and owin object
#wpNGA_ppp = wp_ppp[nga_owin]
wpNGA_ppp_functional = wp_ppp_functional[nga_owin]
wpNGA_ppp_nonfunctional= wp_ppp_nonfunctional[nga_owin]#plot(wpNGA_ppp)
par(mfrow=c(1,2))
plot(wpNGA_ppp_functional)
plot(wpNGA_ppp_nonfunctional)
4.4.2 tmap plots
We can further plot our water points (functional & non-functional) using tmap().
Lets put both functional and non functional water points together. Also, Set the base map to be “OpenStreetMap”.
tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(NGA_wp) +
tm_polygons() +
tm_shape(wp_functional) +
tm_dots(col = "status_clean",
size = 0.01,
border.col = "black",
border.lwd = 0.5) +
tm_shape(wp_nonfunctional) +
tm_dots(col = "status_clean",
size = 0.01,
border.col = "black",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(8, 16))From this point map, even with the different colours points, it is difficult to intrepet the clusters of the functional and non-functional water points. Lets go to the next tab to look at only functional points.
tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(NGA_wp) +
tm_polygons() +
tm_shape(wp_functional) +
tm_dots(col = "status_clean",
size = 0.01,
border.col = "black",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(8, 16))At a macro view (zoom out), there are A LOT of functional water points. While zooming in, it looks like there’s some clusters of water points such as the border between Ife Central and Ife East. Lets take a look at the next tab, wp_nonfunctional.
tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(NGA_wp) +
tm_polygons() +
tm_shape(wp_nonfunctional) +
tm_dots(col = "status_clean",
size = 0.01,
border.col = "black",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(8, 16))Again, at a macro view (zoom out), there are A LOT of functional water points. While zooming in, it looks like there’s some clusters of water points such as the border between Ife Central and Ife East, which seems like the same cluster as shown in wp_functional. This highlights that point map might not suitable in identifying the clusters due to the overlapping points and categorisation. Thus, lets proceed onto the First-Order of Spatial Point Patterns Analysis to plot Kernel density maps which could provide more meaningful insights.
5 First-order Spatial Point Patterns Analysis
Now that we have the point maps from the previous section, moving forward, lets see how would the kernel density maps would display the water points in Osun.
5.1 Kernel Density Estimation
In this section, We will be computing Functional and Non-Functional water points in Osun.
5.1.1 Computing kernel density estimation using automatic bandwidth selection method
We will be using the bw.diggle() method than bw.ppl() because it is more suitable for this exercise to detect single tight cluster in the midst of random noise which was observed earlier section in the point maps, theres seem to be clusters, for example, along the border between Ife Central and Ife East.
kde_wpNGA_functional.bw <- density(wpNGA_ppp_functional,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
kde_wpNGA_nonfunctional.bw <- density(wpNGA_ppp_nonfunctional,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_wpNGA_functional.bw)
plot(kde_wpNGA_nonfunctional.bw)
It looks the unit of measurement is of the default value in meter. That’s why the density values computed is in “number of points per square meter”.
Before moving to the next section, let’s retrieve the bandwidth used to compute the KDE layer.
#bw <- bw.diggle(wpNGA_ppp)
#bw
bw_functional <- bw.diggle(wpNGA_ppp_functional)
bw_functional sigma
252.1687 bw_nonfunctional <- bw.diggle(wpNGA_ppp_nonfunctional)
bw_nonfunctional sigma
308.2061 5.1.2 Rescalling KDE values
wpNGA_ppp_functional.km <- rescale(wpNGA_ppp_functional, 1000, "km")
wpNGA_ppp_nonfunctional.km <- rescale(wpNGA_ppp_nonfunctional, 1000, "km")#kde_wpNGA.bw <- density(wpNGA_ppp.km, sigma=bw.diggle, edge=TRUE, #kernel="gaussian")
#plot(kde_wpNGA.bw)
kde_wpNGA_functional.bw <- density(wpNGA_ppp_functional.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
kde_wpNGA_nonfunctional.bw <- density(wpNGA_ppp_nonfunctional.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_wpNGA_functional.bw)
plot(kde_wpNGA_nonfunctional.bw)
Unlike the point map, from the above KDE maps, it highlight the signifcant clusters and much easier to identify the clusters of functional and non-functional water points clusters.
5.2 Adaptive KDE
5.2.1 Computing KDE by using adaptive bandwidth
Fixed bandwidth method is very sensitive to highly skew distribution of spatial point patterns over geographical units for example urban versus rural. Hence, for this exercise, we will be using adaptive bandwidth instead.
#kde_wpNGA_adaptive <- adaptive.density(wpNGA_ppp.km, method="kernel")
#plot(kde_wpNGA_adaptive)
par(mfrow=c(1,2))
kde_wpNGA_functional_adaptive <- adaptive.density(wpNGA_ppp_functional.km, method="kernel")
plot(kde_wpNGA_functional_adaptive)
kde_wpNGA_nonfunctional_adaptive <- adaptive.density(wpNGA_ppp_nonfunctional.km, method="kernel")
plot(kde_wpNGA_nonfunctional_adaptive)
5.2.2 Converting KDE output into grid object
To make the KDE output is suitable for mapping purposes.
#gridded_kde_wpNGA_bw <- as.SpatialGridDataFrame.im(kde_wpNGA.bw)
#spplot(gridded_kde_wpNGA_bw)
par(mfrow=c(1,2))
gridded_kde_wpNGA_func_bw <- as.SpatialGridDataFrame.im(kde_wpNGA_functional.bw)
spplot(gridded_kde_wpNGA_func_bw)
gridded_kde_wpNGA_nfunc_bw <- as.SpatialGridDataFrame.im(kde_wpNGA_nonfunctional.bw)
spplot(gridded_kde_wpNGA_nfunc_bw)
5.2.2.1 Converting gridded output into raster
Next, we will convert the gridded kernal density objects into RasterLayer object by using raster() of raster package.
#kde_wpNGA_bw_raster <- raster(gridded_kde_wpNGA_bw)
#kde_wpNGA_bw_raster
kde_wpNGA_func_bw_raster <- raster(gridded_kde_wpNGA_func_bw)
kde_wpNGA_func_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045 (x, y)
extent : 176.5032, 291.0438, 331.4347, 454.5201 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : v
values : -5.092436e-15, 25.49435 (min, max)kde_wpNGA_nfunc_bw_raster <- raster(gridded_kde_wpNGA_nfunc_bw)
kde_wpNGA_nfunc_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045 (x, y)
extent : 176.5032, 291.0438, 331.4347, 454.5201 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : v
values : -3.925434e-15, 20.49412 (min, max)WHY? - crs property is NA during the convert
#projection(kde_wpNGA_bw_raster) <- CRS("+init=EPSG:3414")
#kde_wpNGA_bw_raster
projection(kde_wpNGA_func_bw_raster) <- CRS("+init=EPSG:3414")
kde_wpNGA_func_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045 (x, y)
extent : 176.5032, 291.0438, 331.4347, 454.5201 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
source : memory
names : v
values : -5.092436e-15, 25.49435 (min, max)projection(kde_wpNGA_nfunc_bw_raster) <- CRS("+init=EPSG:3414")
kde_wpNGA_nfunc_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045 (x, y)
extent : 176.5032, 291.0438, 331.4347, 454.5201 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
source : memory
names : v
values : -3.925434e-15, 20.49412 (min, max)Lets display the raster in cartographic quality map using tmap package.
#tm_shape(kde_wpNGA_bw_raster) +
# tm_raster("v") +
# tm_layout(legend.position = c("right", "bottom"), frame = FALSE)
par(mfrow=c(1,2))
tm_shape(kde_wpNGA_func_bw_raster) +
tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"), frame = FALSE)tm_shape(kde_wpNGA_nfunc_bw_raster) +
tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"), frame = FALSE)5.2.3 Spatial Point Patterns using KDE
5.2.3.1 Extract Osun
In this section, we will be comparing KDE of water points at Osun as the study area.
#osun = NGA_sc[wp_sc$clean_adm1 == "Osun",]
osun = NGA_sc
plot(osun, main = "Osun")
5.2.3.2 Converting the spatial point data frame into generic sp format
osun_sp = as(osun, "SpatialPolygons")5.2.3.3 Creating owin object
osun_owin = as(osun_sp, "owin")5.2.3.4 Combining wp points and Osun (Study Area)
wp_osun_ppp_func = wp_ppp_functional_jit[osun_owin]
wp_osun_ppp_nfunc = wp_ppp_nonfunctional_jit[osun_owin]Next, rescale() function is used to trasnform the unit of measurement from metre to kilometre.
wp_osun_ppp_func.km = rescale(wp_osun_ppp_func, 1000, "km")
wp_osun_ppp_nfunc.km = rescale(wp_osun_ppp_nfunc, 1000, "km")used to plot these four study areas and the locations of the childcare centres.
par(mfrow=c(1,2))
plot(wp_osun_ppp_func.km, main="Osun")
plot(wp_osun_ppp_nfunc.km, main="Osun")
5.2.3.5 Computing KDE
The code chunk below will be used to compute the KDE of these four planning area. bw.diggle method is used to derive the bandwidth of each
par(mfrow=c(1,2))
plot(density(wp_osun_ppp_func.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Osun Functional")
plot(density(wp_osun_ppp_nfunc.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Osun Non-Functional")
Observations
6 Second-order Spatial Point Patterns
6.1 Analysing Spatial Point Process Using G-Function
In this section, the G function measures the distances between any two events and their respective nearest neighbours. Using spatstat package’s Gest() and envelope() to run a Monte Carlo simulation test, which is used to predict the probability of a range of outcomes when the possibility of random variables is present.
Next, we will be conducting a hypothesis test to confirm the observed spatial patterns above. The hypothesis and test are as follows:
Ho = The distribution of non-functional water points in Osun, Nigeria are randomly distributed.
H1= The distribution of non-functional water points in Osun, Nigeria are NOT randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05.
Monte Carlo test with G-function
G_func_osun.csr <- envelope(wp_osun_ppp_func, Gest, nsim = 39)Generating 39 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39.
Done.plot(G_func_osun.csr)
Conclusion From func_osun.csr, it is observed that at 95% confidence interval, the G(r) far above the G(theo) and the envelope. This implies that the functional water points in Osun, Nigeria are clustered. Thus, we can reject the null hypothesis that the distribution of functional water points are randomly distributed.
6.1.1 Non-Functional
Next, Non-functional
G_nfunc_osun = Gest(wp_osun_ppp_nfunc, correction = "border")
plot(G_nfunc_osun, xlim=c(0,500))
Next, we will be conducting a hypothesis test to confirm the observed spatial patterns above. The hypothesis and test are as follows:
Ho = The distribution of non-functional water points in Osun, Nigeria are randomly distributed.
H1= The distribution of non-functional water points in Osun, Nigeria are NOT randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05.
Monte Carlo test with G-function
G_nfunc_osun.csr <- envelope(wp_osun_ppp_nfunc, Gest, nsim = 39)Generating 39 simulations of CSR ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39.
Done.plot(G_nfunc_osun.csr)
Conclusion From nfunc_osun.csr, it is observed that at 95% confidence interval, the G(r) far above the G(theo) and the envelope. This implies that the non-functional water points in Osun, Nigeria are clustered. Thus, we can reject the null hypothesis that the distribution of non-functional water points are randomly distributed.
:::
7 Analysing Spatial Point Process Using L-Function
Lets also investigate whether the spatial distribution of functional and non-functional water points are independent from each other.
Again, we will be conducting a hypothesis test to confirm the observed spatial patterns above. The hypothesis and test are as follows:
Ho = The distribution of non-functional water points in Osun, Nigeria are randomly distributed.
H1= The distribution of non-functional water points in Osun, Nigeria are NOT randomly distributed.
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05.
Monte Carlo test with L-function
#L_func_osun.csr <- envelope(wp_osun_ppp_func, Lest, nsim = 39, rank = 1, glocal=TRUE)7.0.1 Plot the result from L-function
#plot(L_ck.csr, . - r ~ r, xlab="d", ylab="L(d)-r")Conclusion
8 Acknowledgement
Thank you Prof Kam for our IS415 Geospatial Analytics and Applications course materials & resources